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ANALYSIS AND DESIGN OF JUMP COEFFICIENTS IN DISCRETE 
STOCHASTIC DIFFUSION MODELS 


LINA MEINECKE, STEFAN ENGBLOM, ANDREAS HELLANDER, PER LOTSTEDT* 

Abstract. In computational system biology, the mesoscopic model of reaction-diffusion kinetics 
is described by a continuous time, discrete space Markov process. To simulate diffusion stochastically, 
the jump coefficients are obtained by a discretization of the diffusion equation. Using unstructured 
meshes to represent complicated geometries may lead to negative coefficients when using piecewise 
linear finite elements. Several methods have been proposed to modify the coefficients to enforce the 
non-negativity needed in the stochastic setting. In this paper, we present a method to quantify the 
error introduced by that change. We interpret the modified discretization matrix as the exact finite 
element discretization of a perturbed equation. The forward error, the error between the analytical 
solutions to the original and the perturbed equations, is bounded by the backward error, the error 
between the diffusion of the two equations. We present a backward analysis algorithm to compute the 
diffusion coefficient from a given discretization matrix. The analysis suggests a new way of deriving 
non-negative jump coefficients that minimizes the backward error. The theory is tested in numerical 
experiments indicating that the new method is superior and minimizes also the forward error. 
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1. Introduction. The molecular pathways that regulate cellular function are 
inherently spatial. Cells have a high level of sub-cellular organization, such as a con¬ 
fined nucleus in eukaryotes or membrane bound reaction complexes. Macromolecules 
are transported by passive diffusion or active transport, driven by molecular motors, 
between different areas in the cell in order to arrive at the correct location to perform 
their function. For example, many gene regulatory pathways rely on a cytoplasmic 
component, where a signal is propagated from the cell membrane to the nucleus, 
and a nuclear component, where transcription factors bind to DNA to regulate the 
expression of genes. 

On a macroscopic modeling level, the diffusion equation - a partial differential 
equation (PDE) - is used to describe the time evolution of the concentration of a 
population of molecules undergoing diffusion. This is a valid model if molecules are 
abundant. But in cellular regulatory networks, key proteins such as transcription fac¬ 
tors are only present in low copy numbers and the deterministic PDE model becomes 
inaccurate. Experiments [8, 26, 32, 34, 35, 38, 45] and theory [14, 33] have shown 
the importance of accounting for intrinsic noise when modelling cellular control sys¬ 
tems. Consequently, we need spatial stochastic simulation methods, and diffusion in 
particular is described by a random walk. We can distinguish between two levels of 
accuracy. 

On the mesoscopic level we use a discrete Brownian motion to model the jump 
process of the molecules. The domain is partitioned into compartments or voxels. 
The state of the system is the number of molecules of each species in each voxel. 
Molecules can jump between neighboring voxels and react if they are in the same 
voxel. The probability density function (PDF) for the probability to be in a state 
at a certain time satisfies a master equation. If bimolecular reactions are included 
there is in general no analytical solution for the PDF and a numerical solution is 
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difficult to compute due to the high dimensionality of the state space. Instead, the 
stochastic simulation algorithm (SSA) can be used to generate trajectories of the 
system. It was first developed by Gillespie [16, 17] for reactions independent of space. 
Its efficiency has been improved in [5, 15] and it is extended to space dependency with 
a Cartesian partitioning of the space in [7, 20, 22]. A more accurate description is 
the space-continuous microscopic level where individual molecules are followed along 
their Brownian trajectories. Methods and software for this approach are found in 
[1, 6, 23, 41, 50]. 

In this work, we focus on diffusion at the mesoscopic level. The probability per 
unit of time for a molecule to jump from its voxel to a neighbor is obtained by a 
discretization of the Laplace operator in the diffusion equation on the same mesh. 
A mathematically equivalent interpretation is that they are obtained from a dis¬ 
cretization of the Fokker-Planck equation for Brownian motion. The resulting matrix 
is the generator of a Markov process and all the off-diagonal entries, representing 
transition rates, need to be non-negative. On the macroscopic level of determinis¬ 
tic PDEs, requiring non-negative jump coefficients enforces the discrete maximum 
principle [46, 48]. To represent the complicated geometries present in cells (e.g. mi¬ 
tochondria or convoluted membranes), we work with unstructured meshes, meaning 
triangular or tetrahedral meshes in 2D and 3D. Many interesting cellular processes 
happen on the membranes. Using a vertex centred discretization allows us to couple 
diffusion in the bulk to diffusion on the surface of the domain in a straightforward 
way. If a molecule reaches a boundary node we can use the surface mesh for its 2D 
diffusion while bound to the membrane. 

Piecewise linear finite elements on unstructured meshes are used in [9] to obtain 
the jump propensities. Mature software exists for discretizing PDEs with the finite 
element method (FEM) on a given mesh, e.g. [30]. In 2D, mesh generators are usually 
able to provide high quality meshes [10], but in 3D the mesh quality decreases and 
negative off-diagonal elements often appear in the FEM discretization matrix [3, 24, 
27]. This matrix can then no longer be interpreted as the generator matrix to a Markov 
process and thus provide transition rates. For our application in stochastic simulations 
in systems biology, we need to modify this discretization matrix to guarantee non¬ 
negative jump coefficients. In [9] the negative entries are set to zero and the diagonal 
element is recalculated so that the row sum equals zero. This changes the diffusion 
speed and leads to errors in, for example, the time a signal needs to propagate from the 
nucleus to the cell membrane. To address this, we have in previous work [31] developed 
a method that preserves mean first passage times. Another approach to obtain non¬ 
negative jump coefficients is to use the finite volume method (FVM). But as we will 
see, the vertex centered FVM scheme does not approximate diffusion more accurately 
than a filtered FEM discretization for typical meshes, despite positive coefficients. 
These methods make the stochastic simulation of diffusion mesh dependent, but this 
also holds for the accuracy of the numerical solution of the PDE. FEM and FVM 
coefficients have been modified in [4, 13, 42] to be non-negative but they depend on 
the PDE solution making them unsuitable for stochastic simulation. 

In this paper we analyse the error introduced by modifying the discretization ma¬ 
trix to enforce non-negative jump coefficients. Since the concentration of the species 
simulated by the SSA converges towards the solution of the diffusion equation, [28, 29], 
we quantify the error in this deterministic limit. We use backward analysis to find 
the diffusion equation solved by the new discretization matrix. We study two error 
estimates: the backward error describing the difference in the diffusion in the equa- 
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tions and the forward error describing the error in the solutions to the equations. The 
analysis suggests a new method to obtain a non-negative discretization by minimizing 
the backward error. 

In Section 2 we describe the mesoscopic model and how the jump coefficients 
are obtained for unstructured triangular and tetrahedral meshes. In Section 3, we 
develop theory to bound the forward error by the backward error. An algorithm is 
provided in Section 4 for calculating the backward error and then we show how new 
error minimizing jump coefficients can be computed in Section 5. In the experiments 
in Section 6, we analyse the errors numerically, test our new method, observe that it 
also minimizes the forward error in agreement with the estimates in Section 3, and 
discuss possibilities for a practical implementation. Final conclusions are drawn in 
Section 7. 

Vectors and matrices are written in boldface. A vector u has the components iq 
and the elements of a matrix A are Aij. For vectors and matrices ||u|| p denotes the 
vector norm in £ p and ||A|| p its subordinate matrix norm, and ||u||^ = u T Au for a 
positive definite A. The derivative of a variable u with respect to time t is written 
u t . If u(x), x G S4, varies in space, then ||u||^ p = f n ||u|||Jdfi. 

2. Mesoscopic model of diffusion. To model diffusion in discrete space, the 
domain of interest is partitioned into non-overlapping voxels V&, k = 1,..., N. Each 
voxel Vj has a node or vertex with the coordinates Xj in the interior, see Fig. 2.1. If Vi 
and Vj are neighbors, the vertices x^ and Xj are connected by the edge e^-. Molecules 
in a voxel Vj can diffuse by jumps to a neighboring voxel Vi along the edge e^-. 



Fig. 2.1: The primal triangular mesh (black), defining the edges e^, and the dual 
mesh (blue) defining the voxels V*. 


The state of the system is the discrete number of molecules of each species in each 
voxel. Let jjj be the number of molecules of chemical species Y in Vj. The jump rate 
A ji from Vj to a neighboring Vi needs to satisfy the condition 

(2.1) A ji > 0. 

The total jump rate out of voxel Vj is A j = JV i^j \ ji • The next time for a jump 
from Vj is exponentially distributed with the intensity XjVj. Voxel Vi is chosen as 
the destination with a probability proportional to A ji. After a jump, the number of 
molecules is updated and the time for a new jump is determined. This is the Stochas- 
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tic Simulation Algorithm (SSA) by Gillespie [16] to simulate mesoscopic diffusion of 
molecules between the voxels. 


Algorithm 1 Stochastic Simulation Algorithm for Diffusion [16] 
l: Initialize yk,k = 1,..., TV, in the N voxels at t = 0. 

2: Sample the exponentially distributed time tk with rate XkVk to the first diffusion 
event in all N voxels. 

3 : Let tj be the minimum of all tk . If tj > T then stop, otherwise continue. 

4 : For the jump from Vj, sample a jump to Vi with probability Xji/Xj. 

5: Update t := tj , yi and yj. Sample A U with the rate A iyi and A tj with the rate 
A jyj and recompute U — t + A U and tj=t + A tj. Go to 3. 


We will now show how to derive the propensities Aq from a discretization of the 
diffusion equation. Let y be a vector with entries yi, describing the number of Y 
molecules in voxel Vi. The probability distribution p(y,t) for the distribution of the 
molecules is the solution to the diffusion-master equation 

N N 

( 2 . 2 ) EE •My - t*ij)p (y - - A»j(y)p(y,*), 

i =1 j =1 

where /iqq = — 1 , /iqj = 1 and j±ij,k — 0 for k ^ i,j. Calculating the expected value 
y i of the number of molecules in each voxel i leads to a system of ordinary differential 
equations (ODEs) for the mean concentration iq =Vi/\Vi\ in each voxel Vi 

N IV-I N 

(2.3) u it = TTj-rXji u j — Ui Xij , 

3=1 1 j = l 

see [9]. This can be interpreted as a discretization of the diffusion equation 

(2.4) u t = yA u = V • (yV^), x G D, t > 0, 

n • Viz = 0, x G 9U, t > 0, 

u — i/o, x G U, t = 0, 

with the diffusion coefficient 7 and 7 = 7 I. Thus, the jump rates A j and Xji can be 
computed using discretizations of the diffusion equation on the triangular or tetra¬ 
hedral mesh. The space derivatives in (2.4) are approximated in the voxels by D to 
obtain equations for the unknowns iq 

(2.5) Ut = Du. 

A discretization with FEM using piecewise linear Lagrangean basis and test functions 
yields a mass matrix M and a stiffness matrix S. The diagonal A is obtained after 
mass lumping of M. The diagonal elements are Ajj = \Vj\. Then the system matrix 
in (2.5) is 

(2.6) D = A -1 S. 

Let h be a measure of the mesh size. The solution of (2.5) converges to the solution 
of (2.4) when h —>• 0 and the difference between them is of 0(h 2 ). If the off-diagonal 
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elements Dij in D are non-negative, then these are taken as the jump coefficients Xji 
in the SSA in Algorithm ( 1 ) scaled by the volumes of the voxels |V*| and \Vj\ 


(2.7) 


, _n,M 

\Vj 


IV, 


see [9]. The concentrations Hj/\Vj\ computed by the SSA converge in the limit of 
large numbers of molecules to the concentrations in (2.5) by [28, 29]. 

In 2D, the entry Sij corresponding to edge e^- is 


(2.8) 


Sij = sin (ip + 0)/2 sm(cp) sin( 0 ), 


where ip and 6 are the two angles opposing e^- [49], see Fig. 2.1. If ip + 6 > tt then 
Sij < 0 and we can no longer use it to define a jump propensity. A similar condition 
exists in 3D [49]. Mesh generators in 2D are usually able to construct meshes leading 
to positive Sij [10], but in 3D negative off-diagonal entries often occur [24]. The 
extra requirement in systems biology to have non-negative off-diagonal elements in 
the stiffness matrix is a sufficient but not necessary condition to fulfill the discrete 
maximum principle when solving the PDE (2.4) numerically [46, 48]. 

We now present three different methods of modifying the stiffness matrix S or 
the discretization matrix D containing off-diagonal negative coefficients, so that we 
can interpret them as the generator matrix of the Markov process simulated by the 
SSA. The discretization matrix D is modified to D in [9] such that, if < 0, then 

= 0 and Da = — Xqii Tkj, where rii is the number of edges leaving vertex i. This 
method of calculating the jump coefficients by eliminating the negative contributions 
is denoted here by nnFEM, non-negative FEM. Convergence of the solution to the 
equation with the diffusion operator 7 A is lost but non-negative jump coefficients are 
defined. Solving the system of equations 

(2.9) Uht = A _1 Su/ l = Du/,, 


however, can be viewed as a discrete approximation to a perturbed diffusion equation 


(2.10) u t = V • ( 7 Vu). 

The diffusion matrix 7 belongs to M 2x2 in 2D and M 3x3 in 3D and is symmetric and 
should be positive definite for all x. If 7 is only positive semidefinite it has at least 
one eigenvalue equal to zero which means that there is no diffusion along the direction 
of the corresponding eigenvector. This is unrealistic to happen inside living cells and 
we do not consider this case, although the following analysis can be generalized to the 
positive semidefinite case. 

Another option is to choose a straightforward finite volume method (FVM). If the 
boundary dVj of a voxel Vj consists of rij straight segments ( 2 D) or flat faces (3D) 
dVjiS = 1,... ,rij, of length or area \dVji\ with normal n^ of unit length, then 

(2.11) [ S7 - ( / y\7u) dv = [ n • u ds 'S'' • 7 eji(ui - Uj) 

JVj JdVj \\ e ji\\2 

and the stiffness matrix in (2.9) is Sji = Uji • 7e^|<9V :/ i|/||eji|| 2 . The elements in 
S derived from ( 2 . 11 ) are always non-negative and hence the Xji in (2.7) defined by 
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the FVM are non-negative. The N components of u represent the average value of 
u in the voxels. However, the solution of (2.5) may not converge to the solutions 
of (2.4) when the mesh size h is reduced [44], since the approximation in (2.11) is 
consistent with 7 A u only if the mesh is of Voronoi type [12]. But we can again 
interpret Dij = Sij/\Vi\ as a consistent FEM discretization of the perturbed equation 
(2.10). The scheme in (2.11) is a vertex centered FVM. A cell centered FVM is used 
in [ 21 ] to define the jump coefficients. 

In [31], the jump coefficients are chosen to be close to the FEM coefficients in 

(2.6). If is non-negative, then A ji is as in (2.7). If < 0, then the Xji coefficients 
for voxel j are determined such that the mean first exit time Sj from a vertex j in the 
mesh to the boundary dQ is a solution of the system of linear equations 

(2.12) De = -e, e T = (1,1, ..., 1), 

where = Xji\Vj\/\Vi\. The mean first exit time is the expected time it takes for 
a molecule initially at a position inside H to reach dfl. With these coefficients in 
Algorithm 1 , the average of the simulated first exit times from vertices in the mesh 
agree very well in [31] with the ones computed numerically with a FEM discretization 
of (2.12). This method of computing the jump coefficients is based on the global first 
exit time of the molecules and is denoted by GFET. 

3. Analysis. Previously, we presented three methods that can be regarded as 
modifications of the FEM discretization matrix D into a matrix D with non-negative 
jump coefficients. In this section we view D as a FEM discretization of a certain 
perturbed PDE (3.2) below. To quantify the error introduced by the change in the 
diffusion matrix we therefore aim at bounding the difference between the solutions to 
the PDEs 

(3.1) u t = yA u, 

(3.2) u t = V • (yVft), 

for x £ H, with homogeneous Neumann boundary conditions du/dn = du/dn = 0 
for x £ 3H, and initial data uq = uo at t = 0. Here y(x) is a symmetric, uniformly 
positive definite matrix. The mean first exit time used to define the GFET algorithm 
fulfills Poisson’s equation, [36] 

(3.3) —1 = yAe, 
with the corresponding perturbed equation 

(3.4) -1 = V • (yW) 

and homogeneous Dirichlet boundary condition. 

Let iL 1 (H) be the Hilbert space of all functions u £ for which 

weak derivative exists and lies in L 2 (H). The corresponding weak problems 
and (3.2) are find u,u £ iL 1 (H) such that for \/v £ iL 1 (0), 

(3.5) (v, u t ) = -(Vv, yVu), 

(3.6) (v,ut) = -(Vv,yVu). 

For finite element solutions in a finite dimensional subspace C iL 1 (H) we 

have the set of equations (see also (2.5) and (2.9) above) 

(3.7) Mu t = Su, 

(3.8) Mut = Su. 


the first 
for (3.1) 
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We shall require the following basic a priori estimates 
Lemma 3.1. For some C > 0, 

(3.9) ||Viz || L 2 < \\Vuo\\l 2 , 

(3.10) ||Vix||x ,2 < C\\\7uq\\ L 2, 

where (3.10) assumes that^f is uniformly positive definite such that g ||y ||| < y T y(x)y < 
G||y|||, for some positive constants (g,G), Vx G ft C Vy G where d is the 
dimension. 

Proof The case (3.9) of a scalar 7 follows straightforwardly so we focus on (3.10). 
Letting v = u t in (3.6) we arrive at 

IH ||2 = -(Vw t , 7 Vu) = -“(Vu, 7 Vu), 

since 7 is symmetric. Integrating we get 

( \ 7 u , 7 \ 7 u ) < (Viio,7Viio). 

Invoking the definiteness of 7 we arrive at 

g\\Vu \\ 2 L 2 <G\\Vu 0 \\ 2 L 2 . 

□ 

We now bound the error in the two solutions u and u. 

Theorem 3.2. For some constant C > 0, 

(3.11) ||u - uf L2 < Ct\\qf - 7 ||oo||Vuo||i 2 , 
where the norm of the difference of the diffusion rates is defined by 

(3.12) ll7-7l|oo := max|| 7 - 7 || 2 . 

xEii 


Proof Subtracting (3.5) from (3.6) and using v = u — u we get 

\ f h\\h = -7llVvlll, + (Vv, (7 - 7)V«) < J \ Vt) T (7 - 7)'v«|dfi 

< f Il7 — 7 l| 2 ||Vt;|| 2 ||Vu|| 2 dfl < max || 7 — 7|| 2 f ||Vt;|| 2 ||Vfi|| 2 dfi 

Jn Jn 

< ||7 — 7lloo||Vu||i2 ||Vu||l 2 < ||7-7lloo(||V«|| La + ||Vti|| La )||Vti|| La 
<C||7-7 ||oo||V U o||| 2 

using Lemma 3.1. The estimate (3.11) follows by integration of the inequality. □ 
This shows that the forward error \\u — u\\l 2 is bounded. Using the maximum 
norm of the difference between the two diffusion constants as in (3.11) is, however, 
pessimistic and we now instead use the mean value of ||7 — 7 (x )||2 over U to bound 
the error in the solutions. 

Proposition 3.3. For some constant C > 0, 

(3.13) f\\u - u \\ 2 l2 < C||7 - 7||*(I|Vm||l 4 + ||Vfi|| L 4)||V«|| L <, 
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where 


Il7 - 7II* : = p- Il7 - 7lll dV- 


(3.14) 

Proof. As previously subtracting (3.5) from (3.6) and using v = u — u we get 

\j t M 2 L* < J a 117 - 7l| 2 ||Vt;||2||Vfi|| 2 dfi < (jT || 7 - 7lli^) * ||Vu||l||Vfi||idfi) * 

1 

< Il7-7||*N^ (d f ||v«||l||Vfi||2dnV < ||7-7IW^I^^I|VH|L4||V«|| L 4. 


For t > S > 0 the diffusion equations (3.1) and (3.2) smooth out irregularities in 
the initial data, so that we can assume ||Vix||x ,4 and ||Vft ||£4 are bounded, such that 
for some C > 0 

(3.15) ||Vw || L 4 < C and \\Vu\\ L * < C. 
and hence by integration of (3.13) 

(3.16) || u - u \\ 2 l2 < Ct\\-f - 7 1)*. 

In summary, Theorem 3.2 shows that the forward error can be bounded in terms of 
the difference ||7 — 7 ||oo and some factors that are independent of 7 (assuming that 
7 is uniformly positive definite). Equation (3.16) shows a sharper bound in terms of 
||7 — 7 1|* at the cost of factors that may depend on 7 . On balance we take (3.16) as 
the basis for our further analysis, thus assuming essentially that C in (3.15) depends 
only mildly on 7 . 

The following proposition proves that u and u have identical steady state, which 
shows that the t-dependent estimates in Theorem 3.2 and (3.16) are pessimistic and 
give a relevant bound only for small t. 

Proposition 3.4. Fort — 00 the steady state solutions of (3.1) and (3.2) fulfill 

(3.17) u 00 = u 00 = \\u 0 \\ Ll /\n\. 


Proof Using v = Uoq in (3.6) we have at the steady state that 

0 = (V Uoo 1 7^ Roo) • 

By the positive definiteness of 7 this means Vfioo = 0 and hence Uoq is constant. A 
similar argument for u implies the same property and, moreover, since u is a density 
we can safely assume uo = uo > 0. Setting v = 1 in (3.6) we conclude 

f^ool^l = / Uoodit = (uqo, 1) = (fio? 1) — ll^ollz, 1 
Jn 

and analogously for u oo- □ 

Using uq = uq we have from Proposition 3.4 that, since \\u — u\\l 2 is continuous 
in time, there exists a t* G ( 0 , 00 ) where the error reaches its maximum ||iz(t*) — 
u(t*)\\ L 2 > || u(t) -u(t) || L 2 , Vt. 
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We obtain a similar result for the error in Poisson’s equations. 

Theorem 3.5. Assume dCl G C°° and 'y(x) G C°°, then for the weak solutions e 
and £ of the weak problems corresponding to (3.3) and (3.4) 

(3-18) ||£-£|| 2 L 2 <C|| 7 -7ll* 

for some constant C > 0. 

Proof Choosing v = e — e and subtracting the weak formulations for (3.3) and 
(3.4) we obtain, with ||v||^ = (v, 7 v), 

0 = ||W||? + (Vv, (7 - 7 )Ve), 

using the same arguments as in the proof of Proposition 3.3 gives 

||Vu||? < CH 7 — 7||*I|Vv||l 4 l|Ve||i4. 

By the Poincare-Friedrich’s inequality and the positive-definiteness of 7 
\\v\\l> <C 5 - 1 ||7-7l|.||Vv|| L 4||Ve|| L 4 

By [11, Chap. 6.3, Theorem 6] and the assumptions we obtain 5 G and e G 

This bounds both ||Ve:||oo and ||Vs||oo and hence || < C and ||W || L 4 < 

C. □ 

We conclude that we can effectively bound both the forward error \\u — u\\l 2 and 
the error in the mean first exit time \\e — £\\l 2 by the difference ||7 — 7 ||*- In the 
following section we present an algorithm for calculating this quantity for a given 
discretization matrix D. 

4. Backward analysis. In Section 2, we presented the FVM and the modified 
FEM to compute the stiffness matrix S leading to non-negative jump coefficients in 
(2.7). This matrix can be interpreted as the FEM matrix of a standard, convergent 
discretization of the perturbed equation (3.2). The general diffusion matrix j(x) 
is symmetric, positive definite and may have non-zero off-diagonal elements. The 
difference between 7 and 7 should be as small as possible. This is a measure of how 
close the jump coefficients are to modeling stochastic diffusion which converges to 
isotropic diffusion with a constant 7 . 

4.1. The FEM discretization. Interpreting S as the standard FEM stiffness 
matrix to the perturbed equation (3.2) implies that 

(4.1) Sij = -(V'0 i ,7(x)VV’j) 

for all edges e^-. The sparsity pattern of S and S is the same determined by the 
connectivity of the mesh. Here ^ and fjj are the hat functions of linear Lagrangean 
finite elements with V^( x i) = 1 an d ( x j) = 0 when i 7 ^ j. Since the right hand 
side of (4.1) is a symmetric expression in i and j the perturbed stiffness matrix S has 
to be symmetric. The FVM and nnFEM generate symmetric stiffness matrices. To 
symmetrize S resulting from GFET we use its symmetric part (S + S T )/2 as S in the 
following. The boundary dQ of the domain H is assumed to be polygonal and H is 
discretized such that 


n= (J 

T k er 
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(4.2) 


where T is the set of all non-overlapping elements T k . These elements are triangles 
in 2D and tetrahedra in 3D in the primal mesh on Q defined by the edges e^-, see 
Fig. 2.1. The dual mesh on D defines the voxels Vi in Section 2. With %j C T being 
the set of all triangles in 2D or tetrahedra in 3D containing edge we can write 
(4.1) as 

(4.3) 4 = - Y f VV’f7(x)V^dx=- Y WfL / 

T k eTij Jt « T k eTn Jt * 

= - £ V^ T | Tfc 7.V^| T JT,|, 

since the gradients are constant in T k . It is only the average ~/ k of ^(x) on each 
element T k that contributes to Sij. Thus, we calculate ^ k of the following type in 2D 
and 3D, respectively, 


(4.4) 
With 

(4.5) 


ll D 


f 7/ci 

7/c3\ 

-3D ( 

\%3 

%2j ’ 

7 k = 1 




7fci 

7/c4 

7fc5\ 

7/C4 

%2 

7/e6 1 

7/c5 

7/c6 

77c3 / 


/v^A 

V^2 , 

VVAa ) 


and the coefficients Cijki and L as in Table 4.1, (4.1) becomes for each edge e^- 

L 

(4.6) Sij = ^ y^jCjjkijki. 

T k eTij i=i 



2D 

3D 

L 

3 

6 

k i 
Cijk2 
CijkS 
Cij kA 
Cijk5 
Cijkd 

-VfeV^-2 

-(V^iV^'2 +V^ 2 V^'i) 

-VipnVipji 

-VV’i2VV’j2 

«W> i3 V^3 

-(v^iiV^ +vV’i2VV’ 2 i) 
-(WiiVVh-3 + V^isVV’ii) 
-(V^ 2 VV>j3 + VV’a V^- 2 ) 


Table 4.1: Coefficients for in (4.6). 


In 2D, the integrand in (4.1) is non-zero on two triangles and on at least three 
tetrahedra in 3D for the edges in the interior of Q. One can show using induction 
that a 2D mesh with N vertices and Eb edges at the boundary has T = 2TV — 2 — Eb 
triangles and E = 3 N — 3 — Eb edges. Taking into account that there are three 
unknowns per triangle in (4.4) and one equation (4.6) per edge we have to solve an 
underdetermined system for any triangulation containing more than one triangle with 
3N — 3 — 2 Eb remaining degrees of freedom. 
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In 3D, the system of linear equations defined by (4.6) is also underdetermined if 
the mesh consists of more than one tetrahedron. Each edge in the mesh is an edge of 
at least one tetrahedron but there may be only one tetrahedron associated with the 
edge on the boundary. Then the number of unknowns is six in (4.4) and the number 
of linear constraints (4.6) is six. For each additional tetrahedron sharing the same 
edge, there are six new unknowns and three new constraints. The total number of 
unknowns for each edge is 6 T where T is the number of tetrahedra with a common 
edge and the number of linear constraints is 3 T + 3. Locally, the diffusion matrix 7 
is underdetermined with 3 T — 3 degrees of freedom. 

Consequently, the diffusion 7 satisfying (4.6) is not unique, but for all possible 7 
the error analysis in Section 3 holds. We obtain the sharpest bounds on \\u — u\\l 2 
and \\e — e \\^2 by finding 7 satisfying (4.6) and minimizing the difference H 7 — 7 ||* 
in the equations. An alternative would be to replace || • JI 2 in (3.14) by the Frobenius 
norm || • \\p. Since 

(4.7) ||A|| 2 < ||A|| f < Vd||A|| 2 

for a matrix A [18], the bound in Section 3 is sharper if the minimization is made 
in the || • H 2 norm. In the following, we propose a global and a local optimization 
procedure to find these minimizers 7 . 

4.2. Global optimization. The diffusion matrix 7 closest to the original diffu¬ 
sion 7 with constant coefficient 7 is found by minimizing the distance between 7 (x) 
and 7 under the constraints in (4.6). The stiffness matrix S is given by (2.7) and one 
of the methods in Section 2 . As only the average 7 ^ of 7 (x) appears in the FEM 
approximation on each triangle T&, see (4.3), the norm of the difference in diffusion 
in (3.14) reduces to the weighted sum of the differences \\j k — 7||| as a measure of 
the distance resulting in the following optimization problem 

( 4 - 8 ) min ^ |T fe ||| 7 fc - 7 ||^ 

-A 

T k er 

L 

(d.9) ^ ^ 'y ^ Cjjki'lki = Sij , Ve^f. 

T k ETij 1=1 

This is a nonlinear programming problem with 3 T variables in 2D and 6 T variables 
in 3D and E linear constraints. 

The difference || 7 fe — 71|i * s a convex function in the unknowns in 7 ^. Hence, the 
objective function in (4.8) is a convex function too. Since also the constraint set in 

(4.9) is convex, the local solution to (4.8) and (4.9) is the unique global optimum. 
If Sij > 0 for all i, j from the FEM discretization with diffusion constant 7 , then 
Sij = for all i,j and the solution to (4.8) is 7 ^ = 7 for all T k . 

The mean value matrix 7 ^ defines two (three) main axes in 2D (3D) on T k . Let 
the columns of V be the eigenvectors Vj of 7 fe with eigenvalues A j . After a coordinate 
transformation from x to y with x = Vy the diffusion term is 

(4.10) V x -(jV xU ) = £;A j ~. 

j i 

The eigenvectors define the main axes of the diffusion and the diffusion speed along 
those axes is given by the eigenvalues of 7 . Since 

(4.11) || 7 fc- 7 lb = max|A i - 7 |, 

J 
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the £2 norm in (4.8) measures the maximum deviation in speed of the diffusion in 7 
compared to 7 weighted by the size of T k . In the Frobenius norm 


(4.12) 


Il7 fc -7 ||f 


E(A,-7) 2 j 


and the norm is equal to the £2 norm of the difference in diffusion speed in all direc¬ 
tions. The objective function in (4.8) is continuous in 7 but it is not continuously 
differentiable everywhere. 

4.3. Local optimization. The optimization problem in the previous section 
may be computationally expensive but it is simplified if we approach the solution of 
(4.9) by local optimization. Let Sij be defined by 

(4.13) Sij = {e mn : e mn is an edge of any T k G Tij} 


The adjacent 7 fc in X k in %j for each edge e^- are optimized, while keeping Sij 
constant on the other edges in £ij. Update 7 fc with the most recently computed 
diffusion matrix. Then iterate over all edges once. Still, the underdetermined system 
(4.9) will be satisfied but with a different 7 L compared to 7 ° solving (4.8). The 
algorithm is as follows 


Algorithm 2 Local Optimization I 


1: 7 fe = 7, VT fc G r 

2: for all eij do 

3: Solve 


min^neiy 

^^1^1117^-7111 



^T k eTij 5D;=i Cijki7 k i w 

!CO 

II 


Etl Cmnkatr 

— yVx Cmnkl'lkl 

V &mn ^ Sij \ 6ij , ST k G Tij 

a ^ _ -new rri 

4: Ik ~lk > 1 k 

5: end for 

e Tij 



The diffusion *y k changes successively only on the elements adjacent to e^- (two 
triangles in 2D and at least three tetrahedra in 3D in the interior) in each iterative 
step. At each inner edge in 2 D, there are six variables and five constraints. As 
remarked in Section 4.1 above, the number of variables in 3D in Algorithm 2 is 6 T 
and the number of constraints is 3 T + 3, where T is the number of tetrahedra sharing 
the common edge e^-. For a boundary edge, the number of unknowns equals the 
number of constraints in 2D and one has to solve only the linear system in (4.6). 

Also here we have that if Sij = Sij > 0 for all i,j then the solution is = 7 . 
The order in which the edges are traversed matters for the result by the algorithm but 
when all edges have been visited then (4.9) is satisfied. In the numerical experiments 
in Section 6 , the order is random but other choices are possible. 

The global 7 ^ from (4.8) and the local 7 ^ from Algorithm 2 fulfill 


(4.14) 









since 7 ^ is the global minimum solution. 

Neither the global nor the local procedure to determine 7 guarantee its positive 
definiteness when the solution is computed with an optimization algorithm for a non¬ 
linear objective function with linear constraints. Extra nonlinear constraints can be 
added to enforce positive definiteness. That leads to slow algorithms or sometimes 
very large backward errors H 7 — 7 /JI 2 in the numerical experiments in Section 6 . An 
alternative would be to apply a computationally more expensive semi-definite pro¬ 
gramming algorithm [ 2 , 47] to the problem. In Section 6 , we first compute 7 without 
constraints for positive definiteness and then check the solution for positive definite¬ 
ness. The nonlinear programming algorithm finds positive definite 7 fe for all elements 
in most cases. 

The diffusion 7 is computed for a given mesh of finite mesh size h. What happens 
with 7 when the mesh is refined depends on the mesh generator. If all Sij become 
non-negative as h 0 , then 7 — 7 . Otherwise, there will be a difference ||7 — 7 ||* 
of 0 ( 1 ) as h vanishes. 

The backward analysis is extended in the next section to the design of the stiffness 
matrix S such that the backward error ||7 — 7 ||* is minimized. 

5. Design. In the previous section, we described how to analyze existing meth¬ 
ods for creating positive jump coefficients by backwards analysis. In this section we 
determine a new discretization using FEM by minimizing the backward error. We 
devise non-negative jump coefficients such that the perturbed diffusion 7 is as close 
as possible to the original diffusion with a constant 7. The connectivity of the net¬ 
work of edges is the same as in Section 4.1 but is free to vary. Molecules in the 
stochastic setting are allowed to jump only to the neighboring voxels but the rate is a 
free variable to be optimized such that the distribution of molecules converges to the 
diffusion equation (3.2) in the limit of large molecules numbers. 

5.1. Global optimization. The diffusion 7 ^ in each triangle or tetrahedron is 
determined such that 

(5-1) min J2 \ T k\\\lk-l \\2 

-A 

T k eT 

L 

(5.2) E E Cijkilki > 0, Meij. 

T k (z7ij 1=1 

The equality constraints in (4.9) are replaced by the inequalities in (5.2). The new 
jump coefficients Xji are computed by the optimal 7 and Sij 

L 

(3-3) Sij = ^ ^ ^ ^ Cjjki^fki•> i 7 ^ 3-> 

T k E7\j 1=1 

inserted into (2.7). The stiffness matrix is thus obtained from the FEM discretization 
of the diffusion term in (2.10) with linear Lagrangean elements and diffusion matrix 
Ik on T k- 

5.2. Local optimization. The local optimization algorithm in Section 4.3 to 
analyze given jump coefficients is modified in the same way to generate new coefficients 
instead. For each edge, the adjacent diffusion matrices are computed such that they 
are close to 7 and the non-negativity constraint is satisfied for the edges. Instead 
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of keeping the contribution to the other edges constant we let it vary constrained by 
non-negativity. If e mn is an edge in Eij \ e ^, then we allow 7 to be such that 


(5.4) 


L 


^ ^ ^ ^ Crnnkllfkl — $mn — 0- 

T k eTmn l = 1 


In each local optimization step, in the elements T k adjacent to edge are modified 
while keeping 5 mn in other edges in non-negative. Splitting the sum in (5.4) into 
two parts we have 

L L 

(5-5) £ £ CmnklTkl ^ E E CmnklTkl 

T k eTijnTmn 1 = 1 TkeTmnVTijnTmn) 1 = 1 

L 

(5-b) = ^ ^ ^ ^ CmnklTkl &mn • 

TkeTijnTmn 1 = 1 

The diffusion matrix on the left hand side of (5.5) is updated given the diffusion 
matrix in the right hand side in (5.6). This is repeated successively for all edges in 
the following algorithm. 


Algorithm 3 Local Optimization II 

1: 7 k = 7, VTfc G T 

2 : for all do 

3: = X ]r k eTij ^2i=i CijkHki 

4: end for 
5: for all do 

6 : Solve 

£ Tfcer jT fc | 117^-7111 

Tner^TtiTki^r >0 

"22TkeTijnTmn 22l= 1 C'mnkn 1 k[ W — 22T k eTijnTmn 221=1 CmnklTkl ~ Mmn 

^^mn ^ &ij \ &ij i ^T]k £ 'Tij 

7: lk=Tk eW ,T k £Tij 

8 : Sij = YlT k eTij 22l=1 CijklTkl 

9: end for 


The number of optimization problems to be solved in Algorithm 3 is the number 
of edges E which is bounded by a constant times the number of vertices A in a 
mesh. The size of each optimization problem is independent of E and N . Hence, the 
computational work is proportional to N in the local optimization and of the same 
computational complexity as the matrix assembly of S. The edges are traversed in 
a random order in the experiments in Section 6 . When all edges have been visited 
once, Sij satisfies (5.3) and Sij > 0 and the new Xji is computed using (2.7). 

5.3. Practical implementation. The local minimization problem only con¬ 
tains the adjacent triangles or tetrahedra and is hence faster to compute but rj 2 > 
see (4.14). Instead of running the local Algorithms 2 and 3 only once, we can repeat 
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them iteratively with the results 7 k of the previous iteration as the initial guess for 
the next minimization. Then rj 2 will approach 77 ^?. 

A possibility to speed up the computation is to replace the i 2 norm of the error 
by the Frobenius norm 


(5.7) 


VF 



E iw* 

T k er 


Till- 


The global non-linear minimization problem (4.8) then simplifies to the quadratic 
programming problem 


(5.8) 

(5.9) 


min 7 T H 7 — 2 f T 7 

7 

L 

E E CijkHki> 0, Ve^j, 

T k ETij 1=1 


where H is a diagonal matrix with positive elements on the diagonal and 7 is a vector 
with 7 ki as components. By the relation between the £2 and Frobenius norms (4.7), 
the resulting tjf only yields an upper bound on the global minimum 77 ^. In the local 
optimizations in Algorithms 2 and 3, || • ||| is then substituted by || • |||,. 

We can further reduce the computational complexity by rewriting the high di¬ 
mensional minimization problem (5.8) as the smaller dual problem. 

(5.10) min /x T H/x + 2 f T /x, 


where /x > 0 is equivalent to /x* > 0,Vi, and 

(5.11) H = CH _ 1 C t , f = CH _ 1 f, 7 = -H _ 1 (C t /x - f). 


In (5.11), C is such that (5.9) is replaced by C 7 > 0. The primal problem of dimension 
3 T is hence reduced by approximately a factor two to the dual problem of dimension 
E in 2D. In 3D the dual problem is more than a factor 4.5 smaller than the primal 
problem in numerical experiments in Section 6 . The interior point algorithm is well 
suited for the quadratic programming problems (5.8) and (5.10), see e.g. [2, 25]. 

5.4. Alternatives to determine a non-negative S. Another two possibilities 
are investigated to calculate a S with only non-negative off-diagonal entries. From 
the set of discrete equations (2.5) it appears that a smaller difference between the 
discretization matrices D and D leads to a smaller error in the solution ||u — u\\l 2 - 
That suggests to find a S with the same sparsity pattern as S but with only non¬ 
negative entries such that ||D — £>||2 is minimized. 

A second alternative to guarantee non-negative jump coefficients is adding arti¬ 
ficial viscosity to the system. The same viscosity is added patchwise in all elements 
with a common vertex. If edge e^- corresponds to a negative entry, then enough vis¬ 
cosity to eliminate the negative entry Sij is added to all edges originating from and 
Xj as in the graph Laplacian. The symmetry of the original matrix S is preserved by 
adding | Sij \ /2 to the nodes around x^ and Xj in the following way 

Sik = Sik + \Sij\/2, S ki = S k i + \Sij\/2, Vx/c connected to x* by e fei , 

Sjk — ^jk T \^ij l/^, *5 £j — ^kj H - l^ijl/ 2 , Vx/j connected to x^ by e^. 


This is a generalization of the nnFEM approach where a sufficient amount of viscosity 
is added only to the negative edge. This type of artificial viscocity is introduced in 
[19] to prove that the maximum principle is satisfied for a conservation law. 
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6. Numerical experiments. In this section, we determine numerically the local 
rj L and global rj G backward errors in (4.14) for the different methods generating non¬ 
negative coefficients as described in Sections 2 and 5 with a diffusion coefficient 7 = 1 . 
By the analysis in Section 3, the backward error bounds the forward error of the mean 
values in the spatial distribution of the copy numbers of the molecules \\u— u\\l 2 (3.13) 
and the exit times ||e: — e\\l 2 (3.18). All computations are done in Matlab using its 
optimization routines. The meshes are generated by COMSOL Multiphysics and the 
FEM matrices are assembled by the same software. 

6.1. Diffusion in 2D. The square [—0.5, 0.5] x [—0.5, 0.5] is discretized into 227 
nodes, see Fig. 6.1. As mentioned in Section 2 , mesh generators usually produce good 
quality meshes in 2D and the mesh in Fig. 6.1 is intentionally perturbed to obtain 47 
edges with negative jump coefficients marked by red in Fig. 6.1. 

The requirement to obtain a non-negative discretization poses other constraints 
on the mesh than what is necessary for a FEM solution of high accuracy. Examples of 
quality measures Q related to errors in the finite element solution of Poisson’s equation 
are found in [43]. In 2D, let hi, h 3 , be the lengths of the edges of a triangle of area 
A with the angle (ps opposing the edge of maximum length h 3 . A bound on the error 
in the gradient between / and the approximating fa in the triangle is in [43] 



( 6 . 1 ) 


where Cf is a bound on the second derivatives of /. The measure Q is positive and 
should be as large as possible. Suppose that two triangles with the same edge lengths 
have the edge e^- of length h 3 in common. Then Sij in (2.8) is negative when <^ 3 > tt/2 
while the estimate in ( 6 . 1 ) is as small as possible when <p 3 is in the neighborhood of 
7 r/ 2 . On the other hand, the accuracy is poor if h 3 is large and all angles are less 
than 7 r/ 2 . Then Q is small but the jump coefficients are positive. A large h 3 will of 
course also affect the spatial resolution of the stochastic simulations but the diffusion 
propensities are well defined. 



Fig. 6.1: The mesh in 2D. Negative edges are shown in red. 
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2.5 <e< 3.0 
2.0 <e< 2.0 
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0.6 <e 
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0.3 < e < 0.4 
0.2 <e< 0.3 
0.1 <e< 0.2 
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Local 


Global 


FVM 


GFET 


nnFEM 






MBE 




Fig. 6.2: The backward error calculated by the local and global minimizations in 
Sections 4 and 5. The local error is e = \\~y k — 7 II 2 on each triangle T k G T. 


6.1.1. Backward analysis. We compare the finite volume method (FVM), the 
symmetrized global first exit time method (GFET), the non-negative finite element 
method (nnFEM), and the method minimizing the backward error (MBE). These 
methods all produce discretization matrices approximating the Laplacian with only 
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non-negative off-diagonal entries. The experiments are carried out for the mesh in 
Fig. 6.1. In Fig. 6 . 2 , the local backward error e = \\~y k — 7 H 2 is plotted for all triangles 
Tk G T, calculated by the local and global minimizations in Algorithms 2 and 3 and 
in (4.8), (4.9) and (5.1), (5.2). 

In Table 6.1 we show the error 


(6-2) m= ,/E \Tk\\\ik-i\\im, 

V T k eT 

where 772 is either 77 I' or rj^ according to (4.14). For the global MBE, the matrices 
7 fc are first calculated by nnFEM and then used as an initial guess for the MBE 
optimization. The minima are computed by Matlab’s fmincon with the active-set 
algorithm. 



V 2 

V 2 

FVM 

0.4211 

0.2729 

nnFEM 

0.2057 

0.1524 

GFET 

0.1963 

0.1356 

MBE 

0.0693 

0.0690 


Table 6.1: The global backward errors in (6.2) computed locally 77 ^ an d globally 77 ^. 


The local calculation of the backward error is more pessimistic than the global 
one in the table as expected from (4.14) but the ranking of the different methods 
is the same for both vfe and 77 ^. The FVM naturally leading to non-negative jump 
coefficients causes the largest backward error when used on a poor mesh. A par¬ 
tial explanation to the FVM results may be that the jump coefficients are generated 
by a different principle than FEM. The fluxes over the element boundaries are ap¬ 
proximated in FVM and the method is forced into the framework of FEM. On four 
triangles, a discretization with FVM even leads to a negative definite diffusion matrix 
7 fc when calculated locally without the positive definiteness constraint. The GFET 
and nnFEM perform comparably for the mesh in Fig. 6.1. Computing D for GFET is 
slightly more expensive than setting the negative off-diagonal entries to 0 in nnFEM. 
However, contrary to the nnFEM the GFET preserves the exit time property of the 
original diffusion, see [31]. The minimization constrained by inequalities to obtain the 
discretization matrices with MBE improves the introduced backward error substan¬ 
tially. The faster local and slower global minimization algorithms yield similar results 
for MBE on the mesh in this example. 

6.1.2. Forward error. The relative error in the discrete solution \\uh—Uh\\L 2 /H'M/iIIl 2 
is computed to verify our analysis in Section 3. Here ||^^,H^a = u[M% with being 
the vector of the solution in each node and M the lumped mass matrix. The rela¬ 
tive error is plotted between the discrete solution with the original discretization 
matrix D (with negative off-diagonal entries) and the perturbed discrete solution Uh 
resulting from one of the algorithms generating non-negative off-diagonal entries in D. 

The system of ordinary differential equations in (3.8) is solved by Matlab’s ode 15s. 

Fig. 6.3(a) shows the initial condition 

(6.3) 7i(x, 0) = tanh(20xi) tanh( 20 x 2 ) + 1 
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on the mesh in Fig. 6.1. 

The forward error ||^(x, t) — {^(x, t) ||l 2 /||^/i||l 2 in space at t = 0.01 is depicted 
in Figs. 6.3(b)-(f). 



2 

1.5 





,e 

1.4 


0.8 

0.6 

0.4 

0.2 

0 


(a) Initial Condition 


(b) GFET 


(c) nnFEM 




0.8 

0.6 

0.4 

0.2 

0 


(d) FVM 


(e) Local MBE 


(f) Global MBE 


Fig. 6.3: (a) Initial condition, (b)-(f) Forward error || Uh{x.i) — Uhfai) ||l 2 /IK/i||l 2 for 
different approximations at each node at t = 0.01. 


The forward error behaves in the way predicted by the backward error in Fig. 6.2, 
Table 6.1, and the stability estimates in Section 3. This is also confirmed in Fig. 6.4 
where the forward error in time is displayed in a log-lin scale such that the error for 
short times becomes visible. The unique steady state (3.17) is reached for large t and 
the bound \\uh — Uh\\2 < kt with some k > 0 for the error derived in Section 3 is sharp 
for small t only. 



Fig. 6.4: Forward relative error \\uh — Uh\\ l 2 /\Wh\\ l 2 fc> r 0 < t < 1. 
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Comparing the results in Fig. 6.4 and Table 6.1 we see that the order between the 
methods is the same using the minimization procedure in Section 4 and the solutions 
of (2.5) and (2.9). The performance of the different methods in the forward error 
is correctly predicted by the performance in the backward error as expected from 
Section 3. The MBE is the best and FVM is the worst method but the forward error 
is quite small in all methods with a peak for FVM less than three percent. 


6.1.3. Error in eigenvalues. In the original equation (2.4), the diffusion is 
isotropic and the quotient between the eigenvalues in (4.10) of 7 in 2D is A 1 /A 2 = 1 
and the eigenvectors point in the coordinate directions. The quotient q = A m i n /A max 
is g/G in Lemma 3.1 and for the FVM and the MBE q is found in Fig. 6.5. 
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0 < q < 1/7 
1/7 < q < 2/7 
2/7 < q < 3/7 
3/7 < q < All 
4/7 < q < 5/7 
5/7 < q < 6/7 
6/7 <q < 1.0 
1.0 = q 


(a) 


(b) 


Fig. 6.5: Quotient between the two eigenvalues Xmin/^max of and the direction of 
the eigenvector corresponding to the larger eigenvalue, (a) Global FVM. (b) Global 


MBE. 


Avoiding the negative off-diagonal elements leads to a local anisotropy in the 
diffusion 7 in Fig. 6.5. The eigenvectors corresponding to the larger eigenvalues are 
not aligned but point in what appears to be random directions. The effect of the 
change of the diffusion from 7 to 7 is randomized over D. A global anisotropy is not 
found here in contrast to what we have in the special regular rhombus mesh in [31]. 
All voxels are tilted there in the same direction increasing the diffusion speed in this 
direction in almost all mesh triangles. A random change of the major axis of diffusion 
is expected in general in a mesh created by any mesh generator. 

6.1.4. Alternative methods. The two alternatives suggested in Section 5.4 - 
minimizing the difference between D and D in the £2 norm and adding viscosity - are 
compared to the previous methods. 

At a first glance at (3.7) and (3.8) it may seem like minimizing the difference 
between the matrices D and D in £2 would result in the smallest forward error but this 
is not the case. Indeed, Fig. 6.6 shows the forward error for this approach together with 
the best (MBE) and the worst (FVM) methods. Comparing the matrix deviations 
in Table 6.2 with the results in Fig. 6.6 reveals that there is only a weak correlation 
between a small forward error and the closeness of the matrices in the £2 norm. 
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Viscosity 

FVM 

Local MBE 
Global MBE 
nnFEM 
GFET 
i 2 optimal 

Table 6.2: Relative difference ||D — D|| 


||D-D|| 2 /||D|| 2 

0.7856 

0.4520 

0.1967 

0.1487 

0.1059 

0.0846 

0.0696 

/||D ||2 between the discretization matrices. 



Fig. 6 . 6 : Forward relative error \\uh — Uh\\L 2 /\\ u h\\L 2 for 0 < t < 1. For ^-optimization 
r ]2 = 0.1643. For added viscosity rj^ = 0.6265. 


The I 2 norm of the difference in the discretization matrices in Table 6.2 neither 
reflects the behavior of the forward and backward errors nor does the D optimized in 
the G norm reproduce the correct steady state in Fig. 6 . 6 . There is no unique equation 
and no unique 'y(x) corresponding to a discretization matrix D in Section 4. Hence, 
it is more meaningful to quantify and minimize the error in the solutions Uh and Uh 
than to compare the discretization matrices representing many different analytical 
equations. In Section 3 we showed that the 7 closest to 7 can be used to bound he 
error in the solution. 

We see that adding viscosity to all nodes in the patch leads to adding an unnec¬ 
essarily large amount of viscosity resulting in a larger error than even the FVM in 
this case. 

6.1.5. Practical implementation. The local backward error is too pessimistic 
in Table 6.1, but since only small optimization problems have to be solved involving 
the local triangles or tetrahedra to an edge, it is faster to compute. Furthermore, the 
backward error can be reduced by repeatedly applying Algorithm 2 . The last solution 
in one iteration is the initial guess in the next iteration and the edges are traversed 
in different order in each iteration. 
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Fig. 6.7: The error 77 ^ f° r different number of iterations of Algorithm 2 and nnFEM, 
compared to the error 77 ^ (blue) and the computing time for 772 an d V 2 (green). 


The effect of repeating Algorithm 2 is to spread the local error in each triangle over 
the domain. The 77 ^ error is reduced towards the global minimum, when altering the 
order of the edges in each iteration, but it does not seem to converge to 77 J?. Repeating 
Algorithm 2 about five times in Fig. 6.7 is sufficient to achieve an improved backward 
error at a small increase in computation time. 

Minimizing the £2 norm becomes prohibitively slow especially in the most expen¬ 
sive algorithm of computing the MBE globally. Therefore, to arrive at a practical im¬ 
plementation we switch to the Frobenius norm (5.7). Then the minimization problem 
(5.1) has a quadratic objective function in (5.8) and is solved by Matlab’s quadprog 
with the interior-point-convex algorithm. The approximate computing times in sec¬ 
onds and the computed backward error for the global MBE are displayed in Table 6.3. 


Minimization 

Time 

ijY.neTFkWh-ikWl/M 

Mb 

3620 

0.0690 

IHIf 

0.9737 

0.0804 


Table 6.3: The computation time in seconds for the global MBE and the resulting 
global error in the £2 norm when minimizing the £2 norm and the Frobenius norm in 
the primal (5.8) problem. 


Since || A ||2 < || A||^ the minimization in the Frobenius norm does not reach the 
global minimum in the £2 norm, but this is compensated for by a substantial speedup. 

In Fig. 6.8, we compare the forward errors for the MBE when minimizing in the 
Frobenius and £2 norms. The error resulting from minimization in the Frobenius 
norm is not much larger while a reduction of the computing time of more than 3000 
is achieved in Table 6.3. 
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Fig. 6.8: Forward error of the MBE with minimization in Frobenius and £2 norms. 
The errors for the global and local minimizations in the Frobenius norm are indistin¬ 
guishable. 


6.2. Diffusion in 3D. Our methods are tested on a more realistic mesh such as 
those encountered in systems biology simulations. A sphere with radius 1 is discretized 
into two tetrahedral meshes with 602 and 1660 nodes. In both meshes about 17 percent 
of the edges have a negative jump propensity with the standard FEM discretization. 
Other mesh generators than the one in COMSOL Multiphysics were tested in [24] 
with similar results. 



t 


Fig. 6.9: Forward error \\uh — uh\\l 2 /\\ u h\\L 2 f° r the different methods on a mesh 
with 602 nodes (solid lines) and with 1660 nodes (dashed lines with same markers) 
discretizing the unit sphere. 


In the following experiments, the global backward error and D in MBE are com¬ 
puted by minimizing in the Frobenius norm. The number of unknown variables in the 
global optimization problem for the largest mesh is 51096 in the primal problem (5.8) 
and 10599 in the dual problem (5.10) in Table 6.4. The local optimization problems 
typically have 24 variables in the primal problem and 13 in the dual one. In Fig. 6.9 
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and Table 6.4, we see that the methods examined in 2D behave similarly in 3D. The 
FVM leads to a non-positive definite diffusion 7 ^ in five tetrahedra on the coarse 
mesh and twelve on the fine mesh when calculated locally with Algorithm 2 where 
|| • ||2 is replaced by || • j] f• There is a slight difference between the local and global 
MBE but the error is in general small for all methods. The ranking of the methods 
is the same as in 2 D in Table 6.1 for small t. Since the percentage of negative edges 
is the same when refining the mesh the nnFEM and GFET methods do not improve 
on a finer mesh. How the backward error of the methods behaves when the mesh 
is refined depends not only on the mesh size but also on the shape of the elements. 
For small £, the methods perform in forward error on each mesh as predicted by the 
respective backward error in Table 6.4. 



602 

1660 


Local 

Global 

Local 

Global 

FVM 

0.6415 

0.4284 

0.6077 

0.4227 

nnFEM 

0.3720 

0.2898 

0.3683 

0.2818 

GFET 

0.3604 

0.2618 

0.3570 

0.2632 

MBE 

0.1680 

0.1593 

0.1676 

0.1698 


Table 6.4: Locally and globally computed errors y J2r k eT l^fc|||7fc — Tlli/I^l on th e 
coarse (602 nodes) and fine (1660 nodes) meshes. 


6.3. Mean first hitting time. Molecules in biological cells do not only undergo 
diffusion but also reactions. In order to measure the error in a way relevant for 
reaction-diffusion kinetics, we construct a problem that mimics the mean first binding 
time for two molecules A and B diffusing in a spherical domain in a diffusion limited 
case. The assumption is that the molecules react instantaneously when they are in 
the same voxel. We calculate the mean time it takes for molecule A diffusing in Q 
with reflecting boundary conditions at 9Q to reach a certain node i at where its 
reaction partner B is located. The molecule A is removed when it reaches the node 
at by introducing a sink at this point. This setup models a reaction complex B 
situated at node x^ transforming our molecule of interest A. 

Let pi (x, t ) be the probability distribution function of finding A in x at time t 
when B is at x^ and let S be the Dirac measure. Then pi satisfies 

(6.4) pu (x, t) = 7 Api(x, t) — kS(x — Xj)pj(x, t) 

with a Neumann boundary condition at 90 and a constant k > 0. The mean value of 
the hitting time for A to find B is determined for all possible starting positions of 
A in the mesh. The initial condition is a uniform distribution of A , p^(x, 0) = 1/|0|. 
The domain O is the sphere of radius 1 discretized by the same meshes as in the 
previous section. A discrete approximation of (6.4) is 

(6.5) p it = (7D - Ki)pi, 

where pf = (pn,Pi 2 , • • • ,PiN) and is the zero matrix except for Ku = 10 9 . 

The survival probability Si(t) and the probability density function 7 ^(t) for t* 
are defined by 

Si(t) = [ pi(x,t)dx = P(ji > £), 

Jn 


( 6 . 6 ) 
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The expected value of the hitting time ti can then be calculated by 


(6.7) 


pOO pOO p GO p 

E[r»] = / T7Tj(r) dr = — / rSudr = / / pi(x,T)dxdr 

J o J o J o Jn 

pOC N 

~ / yyVfcl Pikdr. 

fc=l 


In Table 6.5, we compare E[r^] on the coarse mesh for the original discretization 
matrix D and the modified discretizations D described in Sections 2 and 5. A sink is 
placed at one node i in the mesh. Since many interesting reactions in cells occur in 
reaction complexes bound to the membrane or in the nucleus we especially investigate 
the time it takes for A to either find B at a boundary node or at the node closest to 
the center. The average of E[t*] over i is first computed when the sink and B are at 
any node and at any node on the boundary. Then the sink is at the node closest to 
the center. Finally, E[r^] is computed with stochastic simulation employing Algorithm 

l. The time for A to reach B at the center is recorded for each trajectory and the 
average is taken over 10 5 realizations. The inital position of A is sampled from a 
uniform distribution. With N being the total number of nodes, Nb the number of 
boundary nodes, M the number of trajectories, and r ™ the hitting time of trajectory 

m, the quantities in the table are 


1 N x 

Baw = ^E[ri], E'Bnd = ^ E[ri], A’cdet = E[r c ], 


( 6 . 8 ) 


N ^ i onu N 

i= 1 Xi£dQ 

M 


E’Cstoch — M E 


T c 5 ^Std — 


m= 1 


N b 


^ Yi ( E N] - ^Bnd) 2 . 




The results with standard FEM are found in the top row in the table as reference 
values. The FEM values are second order accurate and converge to the analytical 
values of the original diffusion equation when the mesh size is reduced. Since the 
FEM stiffness matrix has negative off-diagonal elements, stochastic simulation with 
its jump coefficients is impossible. In Fig. 6.10, we illustrate the results in Table 6.5 
by plotting the expected time to reach a node as a function of its radial position. 
We average the expected exit times at all nodes in shells of the sphere of width 0.1. 



Eau 

^Bnd 

E’cdet 

E'Cstoch 

FEM 

8.0413 

11.6860 

4.9211 

N/A 

FVM 

8.8255 

12.7055 

5.9304 

5.9733 

GFET 

7.7308 

11.4102 

3.6250 

3.6424 

nnFEM 

7.4794 

10.8482 

4.4648 

4.4707 

MBE 

8.2944 

12.0863 

5.3001 

5.3323 


Table 6.5: Averages of the expected first hitting time E[r*] defined in (6.8) for different 
methods in Columns 2-5 on the mesh with 602 nodes. 
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Fig. 6.10: Expected exit times as function of the radial position. 


The MBE is the superior method both in the average over all nodes and to reach 
the center when compared to the FEM in Table 6.5 and Fig. 6.10. The GFET which 
was designed in [31] to be accurate for the global first exit time - meaning the first 
hitting time of a boundary node - performs best for the boundary nodes. The GFET 
is, however, not able to compute the time to reach the center of the cell very accurately. 
This shortcoming of the method was discussed in [31]. The FVM is too slow with a 
longer time to reach the sinks than FEM and GFET and nnFEM are a little too fast 
corresponding to a global diffusion coefficient larger than 7. This tendency was noted 
also in [31]. The results of the stochastic simulations in Column 5 for the central node 
are close to the deterministic values of (6.7) in Column 4 as expected for a large M. 

The standard deviation of the mean time it takes to reach a boundary node 
measures the variation between different nodes at the boundary and should be small 
(ideally 0) for a good mesh and an accurate discretization. In Table 6.6, we compare 
the expected time to reach a boundary node and its standard deviation (6.8) for the 
different methods on the two meshes. 



602 

1660 


^Bnd 

^Std/^Bnd 

^Bnd 

^Std/^Bnd 

FEM 

11.6860 

0.0876 

15.8015 

0.0808 

FVM 

12.7055 

0.1557 

16.7264 

0.1386 

GFET 

11.4102 

0.1170 

15.4179 

0.1060 

nnFEM 

10.8482 

0.1171 

14.6743 

0.1060 

MBE 

12.0863 

0.0908 

16.2679 

0.0853 


Table 6.6: Expected time to reach a boundary node and its standard deviation on a 
mesh with 602 nodes (left) and a fine mesh with 1660 nodes (right). 


The standard deviation of MBE is close to that of FEM demonstrating that the 
anisotropy in the diffusion introduced by MBE has a small impact compared to the 
accuracy effects of the discretization and the mesh. The small difference in Estd/^Bnd 
between MBE and FEM is most likely explained by the random directions of maximum 
and minimum diffusion as in Fig. 6.5. The relative standard deviation is reduced 
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slightly when the mesh is refined but Eb nd has not yet converged. The Ebnd closest 
to the FEM value is obtained by GFET. When simulating a signal being transmitted 
inside the cell it is advantageous to use GFET if the important reactions occur on 
the membrane. If the signal on the other hand is travelling inside the cytoplasm 
and reacting there, then the MBE results in the most accurate transmission time in 
Table 6.5. 

7. Conclusion. For the discrete stochastic simulation of diffusion in systems bi¬ 
ology, we need jump propensities for the molecules in the discrete space model. These 
propensities are chosen as the off-diagonal elements of the discretization matrix ob¬ 
tained by a numerical approximation of the Laplacian. The jump coefficients have to 
be non-negative. For unstructured meshes, non-negative off-diagonal elements cannot 
be guaranteed with a discretization matrix assembled by a standard finite element 
method (FEM) but there exist different approaches to change this discretization ma¬ 
trix to fulfill the non-negativity condition. As a result of this change, a diffusion 
equation with an altered diffusion is approximated. 

We first present a method to analyze these existing methods producing non¬ 
negative jump propensities on an unstructured mesh of poor quality. The difference 
between the solution to the original and the perturbed diffusion equations is bounded 
by the difference in the diffusion coefficients. Then the perturbed diffusion is retrieved 
by backward analysis. This leads us to the derivation of a new algorithm creating a 
discretization matrix based on FEM, minimizing the backward error. 

We show in numerical experiments that the finite volume method (FVM) to com¬ 
pute a non-negative discretization incurs high forward and backward errors on our 
meshes. Our previously proposed methods of eliminating the negative entries in the 
finite element matrix (nnFEM) and satisfying the global first exit time constraint 
(GFET) perform comparably. The new method to generate jump coefficients on a 
given mesh proposed in this paper results in a considerably smaller error on both an 
artificial mesh in 2D and a realistic mesh in 3D. 

The average of the first hitting time obtained by stochastic simulations with non¬ 
negative jump coefficients is close to the solution of a deterministic equation with a 
modified diffusion as expected. The accuracy of this average compared to the exact 
analytical values not only depends on the number of trajectories in the Monte Carlo 
simulation but also on the mesh size and the mesh quality. In general, with the MBE 
the backward and forward errors are small and the mean hitting time to any node 
is well approximated. The errors in the stochastic diffusion simulation are of the 
same order as the errors in biological measurements [37, 39]. Furthermore, there is a 
variation in the diffusion constant across the cell in measurements in [40] comparable 
to the variation of the space dependent 7 of the perturbed diffusion equation (3.2) 
determined by our algorithms. 

Since the off-diagonal elements are non-negative, the FEM discretization of the 
equation with the modified diffusion satisfies the sufficient conditions for the discrete 
maximum principle for the FEM solution. The discrete maximum principle being 
satisfied for the original diffusion on any mesh, seems to be possible only with a non¬ 
linear scheme as in [4]. Then the stiffness matrix S is reassembled in every time step. 
Our scheme is linear with a constant S but for a modified diffusion. 
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